home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
CD ROM Paradise Collection 4
/
CD ROM Paradise Collection 4 1995 Nov.iso
/
science
/
fastmap.zip
/
HOMOGC.C
< prev
next >
Wrap
C/C++ Source or Header
|
1991-12-09
|
7KB
|
245 lines
/***************************************
* *
* HOMOG *
* *
* Copyright Dr. Jurg Ott 1985 *
* *
* transcribed from Fortran to C *
* by Dr. David Curtis 1988 *
***************************************/
/*
In:
Jurg Ott 1985 Analysis of Human Genetic Linkage, John Hopkins University Press
To carry out the A-test (Ott,Ann Hum Genet 47:311-320,1983) estimating the
proportion alpha of "linked" families jointly with the recombination fraction
theta in these families. Input (free format) and output are interactive or on
files.
Input information
NR = number of lod scores per family (do not include z=0 at rec. freq. 0.5)
ITAB = 1 if lod table is desired, = 0 for test and estimates only
NR rec. freq. (theta) values in increasing order, e.g. .01 .05 ...
Then input the following information for each family:
ISW = 2 if number of recombinants and nonrecombinants are given for next family
= 1 if lod scores follow for next family
= 0 to finish entering family information
V[i] or RK1 RK2 = lod scores or number of recombinants for family. Note that
lods smaller than -30, e.g. -35, are taken to represent minus infinity and
will appear on output as -99.
*/
#include <stdio.h>
#include <math.h>
FILE *fi,*fo;
#define MAXL 20 /* maximum number of lod scores */
#define CONS .30103
#define DLIM -30.
#define ELIM 1.0E-30
#define MAXF 50 /* maximum number of families */
#define STEP .05
int ISW,ITAB,INF[MAXL],NF,i,j,l,NR,NR1;
double RL[MAXF][MAXL],V[MAXL],VS[MAXL+1],THET[MAXL+1];
double AL,R;
double pnorm();
void wr_VS(),comp_res(),get_data();
FILE *get_file();
double pnorm(x2)
double x2;
{
double t;
if (x2<=0) return .5;
t=1/(1+.33267*sqrt(x2));
return .39894228*exp(-.5*x2)*t*(.4361836+t*(.937298*t-.1201676));
}
void wr_VS()
{
int l;
fprintf(fo," %5.2lf ",AL);
for (l=1;l<=NR1;++l)
{
fprintf(fo,"%10.4lf",VS[l]);
if (l==10) fprintf(fo,"\n ");
}
fprintf(fo,"\n");
}
FILE *get_file(mode)
char *mode;
{
char line[240],s[20];
FILE *fp;
do {
*s=*line='\0';
printf("\nEnter %sput file name (just [ENTER] for keyboard) ---> ",mode[0]=='r'?"in":"out") ;
gets(line);
if (sscanf(line," %s",s)!=1) return mode[0]=='r'?stdin:stdout;
} while ((fp=fopen(s,mode))==0 && (printf("\nCould not open %s\n",s),1));
return fp;
}
main()
{
fi=get_file("r");
fo=get_file("w");
fprintf(fo,"\n HOMOGENEITY TESTS\n");
for (;;)
{
get_data();
comp_res();
}
return 0;
}
void get_data()
{
double RKK,RK1,RK2;
printf("\n Enter number of lod scores per family (0 to exit program) ---> ");
fscanf(fi," %d",&NR);
if (!NR) exit (0);
if (NR>MAXL) {fprintf(stderr,"Number of lod scores too large\n");exit(1);}
NR1=NR+1;
printf("\nLod tables desired? (1=Yes, 0=No) ---> ");
fscanf(fi," %d",&ITAB);
printf("\nEnter %d rec. fr. values:\n",NR);
for (i=1;i<=NR;++i)
{
fscanf(fi," %lf",&THET[i]);
#if 0
if (THET[i]>1)
{printf(stderr,"Error: Rec. fr. values must be less than 1\n");exit(3);}
/* Allow user to enter map distances or whatever instead of thetas
- A-test should be just as valid, provided lod scores are entered
and not number of recombinant/non-recombinants.
*/
#endif
INF[i]=VS[i]=0;
}
THET[NR1]=0.5;
VS[NR1]=0;
for (NF=0;;)
{
puts("\nNew family? (2=yes, with no. of rec. and nonrec.;");
printf(" 1=yes, with lod scores; 0=no, carry out test) ---> ");
fscanf(fi," %d",&ISW);
if (!ISW) return;
++NF;
if (NF>MAXF) {fprintf(stderr,"\nNumber of families too large!\n");exit(2);}
if (ISW!=1)
{
printf("\nEnter number of recombinants and number of non-rec. ---> ");
fscanf(fi," %lf %lf",&RK1,&RK2);
RKK=CONS*(RK1+RK2);
for (i=1;i<=NR;++i)
{
R=THET[i];
V[i]=(RK1&&!R)?-90:RKK+RK1*log10(R)+RK2*log10(1-R);
}
}
else
{
printf("\nEnter %d lod scores ---> ",NR);
for (i=1;i<=NR;++i)
fscanf(fi,"%lf",&V[i]);
}
for (i=1;i<=NR;++i)
{
VS[i]+=(R=V[i]);
if (R<DLIM)
{RL[NF][i]=0;INF[i]=1;}
else RL[NF][i]=pow(10.,R);
}
}
}
void comp_res()
{
int JM=NR1,JT;
double VM=0,VT,SUM;
double ALT,CC1,C2,CCC2,E1,E2,E3,AL1;
for (i=1;i<=NR;++i)
{
if (INF[i]) VS[i]=-99;
else if (VS[i]>VM)
{
VM=VS[i];
JM=i;
}
}
CC1=4.60517*VM;
E2=pnorm(CC1);
JT=JM;
ALT=1;
VT=VM;
AL=1;
if (ITAB)
{
fprintf(fo,"\n Lod table\n\n Alpha\n");
wr_VS();
}
while ((AL-=STEP)>.01)
{
AL1=1-AL;
for (l=1;l<=NR;++l)
{
SUM=1;
for (i=1;i<=NF;++i) SUM*=AL*RL[i][l]+AL1;
/* Note: in original Fortran listing */
/* this line should read: */
/* SUM=SUM*(AL*RL(I,L)+AL1) */
/* and not as printed */
if (SUM<ELIM) SUM=-99;
else
{
SUM=log10(SUM);
if (SUM>VT)
{
VT=SUM;
ALT=AL;
JT=l;
}
}
VS[l]=SUM;
}
if (ITAB) wr_VS();
}
if (ITAB)
{
AL=0;
for(l=1;l<=NR;++l) VS[l]=0;
wr_VS();
fprintf(fo,"\n ");for(j=1;j<=NR1;++j) fprintf(fo,"%10.4lf",THET[j]);
fprintf(fo,"\n Theta\n");
}
if (JT==NR1) C2=CCC2=ALT=0;
else
{
C2=4.60517*(VT-VM);
CCC2=4.60517*VT;
}
E1=pnorm(C2);
E3=0.5*exp(-.5*CCC2);
fprintf(fo,"\n Hypotheses Max. lod Alpha Theta\n");
fprintf(fo," H2%18.3lf%7.2lf%7.2lf\n",VT,ALT,THET[JT]);
fprintf(fo," H1%18.3lf (1)%7.2lf\n H0 (0) (0) (0.5)",VM,THET[JM]);
fprintf(fo,"\n Components of chi-square\n\n Source ");
fprintf(fo,"d.f. chi-square p-value\n");
fprintf(fo," H2 vs. H1 Heterogeneity 1%13.3lf%9.4lf\n",C2,E1);
fprintf(fo," H1 vs. H0 Linkage 1%13.3lf%9.4lf\n",CC1,E2);
fprintf(fo," H2 vs. H0 Total 2%13.3lf%9.4lf\n",CCC2,E3);
fprintf(fo,"\n\n Family nr. Posterior prob. of linkage\n");
AL=1-ALT;
for (i=1;i<=NF;++i)
{
R=ALT*RL[i][JT];
R/=(R+AL);
fprintf(fo," %8d%22.4lf\n",i,R);
}
}